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Abstract 

The Ensemble Kalman filter and Ensemble square root filters are data assimila¬ 
tion methods used to combine high dimensional nonlinear models with observed data. 
These methods have proved to be indispensable tools in science and engineering as 
they allow computationally cheap, low dimensional ensemble state approximation for 
extremely high dimensional turbulent forecast models. Erom a theoretical perspective, 
these methods are poorly understood, with the exception of a recently established but 
still incomplete nonlinear stability theory. Moreover, recent numerical and theoretical 
studies of catastrophic filter divergence have indicated that stability is a genuine math¬ 
ematical concern and can not be taken for granted in implementation. In this article 
we propose a simple modification of ensemble based methods which resolves these sta¬ 
bility issues entirely. The method involves a new type of adaptive covariance inflation, 
which comes with minimal additional cost. We develop a complete nonlinear stability 
theory for the adaptive method, yielding Lyapunov functions and geometric ergodicity 
under weak assumptions. We present numerical evidence which suggests the adaptive 
methods have improved accuracy over standard methods and completely eliminate 
catastrophic filter divergence. This enhanced stability allows for the use of extremely 
cheap, unstable forecast integrators, which would otherwise lead to widespread filter 
malfunction. 


1 Introduction 

With the growing importance of accurate weather forecasting and expanding avail¬ 
ability of geophysical measurement, data assimilation for high dimensional dynamical 
system and data has never been more crucial. The ensemble Kalman filter (EnKE) 
[T] and ensemble sqnare root filters (ESRF) O [3] are ensemble based algorithms well 
designed for this purpose. They quantify the uncertainty of an nnderlying system using 
the sample information of a moderate size ensemble thereby significantly 

reducing the computational cost. The simplicity of these algorithms and their accnrate 
performance has fueled their wide application in various fields of geophysical science 

mM- 

Despite their snccess, there are many unresolved issues with ensemble based meth¬ 
ods. First, there is very little theoretical understanding of the methods and notably 
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the stability framework is incomplete. Only very recently has the stability theory 
been partially understood in the finite ensemble size scenario, with progress on well- 
posedness in the fully observed case [6] and subsequently in nonlinear stability of the 
partially observed case, but under the so-called observable energy criterion [7j. A better 
understanding of filter stability is sorely lacking, with recent numerical studies 
revealing the mechanistically mysterious phenomenon known as catstrophic filter di¬ 
vergence, whereby filter state estimates tend to machine infinity, whilst the underlying 
signal remains in a bounded set. In m it has been established rigorously, through an 
elementary model, that this divergence is not caused by the instability of the numerical 
integrator alone, instead the update step in the filter itself plays a crucial role in the 
genesis of divergence. 

In this article we propose a simple modification of EnKF (and ESRF) which resolves 
these issues completely. The modification is a type of covariance inflation, a widely used 
strategy for stabilizing and improving the accuracy of filters. Typically the forecast 
covariance Cn is (in the case of EnKF) additively inflated to Cn + for some choice 
of inflation constant A^- Since the forecast covariance decides how much uncertainty is 
held by the forecast prediction, inflation has the affect of pulling the hlter back towards 
the observations, yielding improved stability. Existing methods of covariance inflation, 
such as constant additive inflation (An = constant) tend to improve accuracy, but are 
still vulnerable to stability issue like catastrophic filter divergence [8]. 

The modification we propose selects the inflation strength A„ adaptively and varies 
according to the distribution of the ensemble. In particular, if the filter is deemed to 
be performing well, the inflation strength is set to zero and the method is reduced to 
EnKF (or ESRF, as desired). If the filter is deemed to be ‘malfunctioning’, then the 
adaptive inflation is triggered. The strength of the inflation becomes larger when the 
hlter strays further into malfunction. To decide when and to what extent the hlter is 
malfunctioning, we employ two simple statistics of the ensemble, 0^ and which 
are based respectively on the ensemble innovation and the cross correlation between 
observed and unobserved components. The two statistics are so chosen as it is clear 
from the theoretical framework that these are precisely the two variables which must be 
controlled to guarantee stability. Nevertheless, there is a quite natural interpretation 
as to why these two statistics are an effective gauge of hlter performance. The full 
derivation and explanation of the adaptively inhated methods are given in Section [3l 

In Sections [Hand El we develop a complete stability theory for the adaptively inhated 
methods by extending the stability framework established in [7j. This framework is 
comprised of two main results: time uniform mean-square estimates on the ensemble 
via a Lyapunov argument and geometric ergodicity of the signal-ensemble process. We 
prove that if the underlying model satishes energy dissipation, then the hlter inherits 
this dissipation and in particular has a Lyapnuov function with compact sub-level 
sets. This is a vast improvement on the results in [7j for EnKF, since hrstly the 
observable energy criterion is no longer required and secondly the Lyapunov function 
for the hlter is guaranteed to have compact sub-level sets. This latter fact leads to 
geometric ergodicity of the signal-ensemble process for the adaptively inhated EnKF, 
which follows as an immediate corollary of the results in [7]. 

In Section [7] we investigate the performance of the adaptively inhated hlter numer¬ 
ically, comparing performance with standard methods, such as non-inhated EnKF and 
non-adaptively (such as constantly) inhated EnKF. The adaptive method performs at 
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least as well and typically better than standard methods and as expected completely 
avoids the issue of catastrophic filter divergence , which is quite prevalent in all stan¬ 
dard methods. The additional expense in computing the inflation strength is minimal 
and hence the adaptive methods run at almost identical speed to the standard meth¬ 
ods. Most impressively, the adaptive method allows for the use of extremely cheap but 
unstable integrators, like explicit Euler. We will see that such integrators are useless 
for standard methods due to the prevalence of catastrophic filter divergence, but are 
very successful when used with adaptive methods. 

Regarding existing literature, there have recently been many methods proposed 
which employ a type of adaptive inflation for ensemble based methods dn Ea [la 
EKE], but our method is truly novel. Moreover, none of the cited methods have an 
established theoretical understanding, and our framework may serve as a good starting 
point for developing such an understanding. 

The structure of the article is as follows. In Section [2 we define the standard en¬ 
semble based methods, EnKF and ESREs. In Section [3] we define the adaptive inflated 
modifications of EnKE and ESREs. In Section H| we derive time uniform mean-square 
estimates for the adaptively inflated filters via a Lyapunov function argument. In Sec¬ 
tion [5] we prove geometric ergodicity of the signal-ensemble process for the adaptively 
inflated filters. In Section [7| we investigate the performance of these filters numerically. 
In Section [8] we conclude with a discussion. The Appendix contains technical tools 
that will be used for the main results. 


2 Ensemble based methods and adaptive infla¬ 
tion 


2.1 Model setup 

In this paper, we assume the signal (truth) sequence Un G is generated by 


Un — ^h{Un-l) + Cn ) 


( 2 . 1 ) 


where T/j : —>• is a deterministic mapping and {Cn}n>i is the system noise. 

The noise Cn is assumed to be is independent of Ci; • • • Cn-i when conditioned on the 
realization of Un-i, the conditional mean is zero E(Cn|L^n-i) = 0 and the conditional 
covariance E(Cn ® Cn\Un-i) is denoted by Rh{Un-i)- Due to the conditional indepen¬ 
dence of the noise sequence, this gives rise to a Markov chain. 

For example, the model is often generated by the solution of a stochastic differential 
equation (SDE) 

dut = 'il^{ut)dt + TidWt , ( 2 . 2 ) 

for a sufficiently regular vector field : R'^ —>• R'^, diffusion coefficient S G and 

e-dimensional Wiener process W. We then take Un = Unh for some fixed h > 0. In the 
notation above, we have 'I'/i(no) = E(u/i|no) and Cn = Unh — hiu(n-i)h)■ It is easy to 
see that this satisfies the above conditions. 
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2.2 Linear observations 


We assume that the truth Un is observed linearly with a mean zero Gaussian pertur¬ 
bation 

Zn = HUn + Cn • 

where the observation matrix H is of size q x d with H = iHQ,0qx(d-q)) where Hq = 
diag(/ii,..., hq), q < d is the rank of H and hi,..., hg > 0 are fixed scalars. For the 
observational noise, we assume that = 0 and (g) ^n\Un-i) = Iq- 

This seemingly restrictive setting can be assumed without loss of generality. Indeed, 
any observational noise covariance can be reduced to the identity via a simple rotation 
on the filtering problem. Suppose that has a nonsingular covariance matrix F (we do 
not consider the singular case in this article) and has an SVD decomposition 

= <I>A'I'^, then we rotate the coordinate system and consider 

Un = ^^Un, L = = Wn + In- (2.3) 

Hence this change of coordinates reduces the observation matrix and the observational 
noise covariance to the desired form. 

If the observation dimension q is larger than the model dimension d, the last d — q 
diagonal entries of A are zero, so the last d — q rows of are independent of the signal 
and play no role in filtering, hence we can ignore and set d = q. 

Since all the transformations above are linear and bijective, filtering Un with is 
equivalent to filtering Un with in the sense that the subsequent assumptions and 
results hold equally for the original and transformed system. When necessary, this will 
be clarified within the text. 

2.3 Ensemble Kalman filter 

In the standard Kalman filtering theory, the conditional distribution of the signal pro¬ 
cess Un given the observation sequence Zi,..., Zn is given by a Gaussian distribution. 
EnKF inherits this idea by using a group of ensembles to represent this 

Gaussian distribution, as the mean and covariance can be taken as the ensemble mean 
and covariance. The EnKF operates very much like a Kalman filter, except its forecast 
step requires a Monte Garlo simulation due to the nonlinearity of the system. In detail, 
the EnKF is an iteration of following two steps, with (for instance) being sampled 
from the equilibrium measure of Un- 

• Forecast step: from the posterior ensemble at time n — 1, a forecast 

ensemble for time n is generated by 

KW = T,(Ki!\) + CW, (2.4) 

(k) 

where Q are independent samples drawn from the same distribution as Cn- 

• Analysis step: upon receiving the new observation Zn, random perturbations of 
it are generated by adding ^n'^: 

v{k) _ fik) 
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(k) 

where ^ are independent samples drawn from the same distribution as Using 
the Kalman update rule, each ensemble member is then updated as follow with 
Cn being the sample covariance of the forecast ensemble: 

uw = uw - - zW) 

= (/ + + (/ + CnH^Hy^dnH^zi^^ , 

where 

A:=l k=l 

See [6] for a derivation of EnKF as an approximate Monte Carlo method for sampling 
the true posterior distribution. In the case of linear model dynamics, it has been 
rigorously shown [MllITl that EnKF is an approximate Monte Carlo method for the 
true posterior distribution, with convergence in the K ^ oo limit. 

Based on our description above, the augmented process {f7„, ..., is 

a Markov chain. As above, we will denote the natural filtration up to time n as 
= <7{Um, ..., m < Ti}, and denote the conditional expectation with 

respect to Tn as E„. 

In the sequel it will be convenient to split the ensemble into observed and unobserved 
components. In particular, for each ensemble member we write Vn^'^ = 
where G and € kei{H). Similarly, we write = {Xn\Yn^'^) and for 
the forecast ensemble covariance 


(2.5) 


( 2 . 6 ) 


a, = 


Bn 

B^, Cn 


Y 


where 


= 


- Xn) ® (X(") - XO ^ 

k=l k=l 


and so forth the other components. In this notation, the update rule ()2.5I1 becomes 
= aW - + Hod^Hl)-\HoX^^^ - Z^)) (2.7) 

for the observed components and 

yP = yW - {BnfH'^{I + HoC^Hl)-\H,xlt^ - ZW) (2.8) 

for the unobserved components. 


2.4 Ensemble square root filters 

(k) 

One drawback of EnKF comes from its usage of artificial noise Cn ^ as this introduces 
unnecessary sampling errors, particularly when the ensemble size is small [18]. The 
motivation behind the artificial noise is to make the posterior ensemble covariance 

C'n := ® 

k=l k=l 
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(2.9) 


satisfy the covariance update of the standard Kalman filter 

Cn = Cn- dnH^iH^CnH + ly^Hdn, 

when the left hand is averaged over da 1201 ED . ESRFs, including the ensemble 
transform Kalman filter (ETKF) and the ensemble adjustment Kalman filter (EAKF), 
aim to resolve this issue by manipulating the posterior spreads to ensure that (I2.9ji 
holds. Both ETKF and EAKF algorithms are described by the following update steps, 
with the only difference occurring in the assimilation step for the spread. As with 
EnKF, the initial ensemble is (for instance) sampled from the equilibrium 

distribution of C/„. 

• Forecast step: identical to EnKF, the forecast ensembles at time n is generated 
from posterior ensembles at time n — 1: 

The forecast ensemble covariance Cn is then computed using (j2.6p . 

• Assimilation step for the mean: upon receiving the new observation the 
posterior ensemble mean is updated through 

_ — - - 1 ^ 

Vn = Vn-CnH^{I + HCnH^)-\HVn-Zn), Vn = (2.10) 

k=l 

• Assimilation step for the spread: The forecast ensemble spread is given by the 

d X K matrix _ _ 

Sn = [VP-Vn,...,V^^^-Vn]. 

To update the posterior spread, first find a matrix E (for ETKF) or 

An E (for EAKF) such that 

-r^TnSn ® TnSn = ® ^ " CnH^{H^CnH + ly^Rdn • 

A — i A — i 

( 2 - 11 ) 

The posterior spread is updated to Sn = TnSn (for ETKF) or Sn = SnAn 
(EAKF), and the ensemble members are updated to 

1/W=F„ + 5W, 


(k) 

where Sn ' denotes the /c-th column of the updated spread matrix Sn- By con¬ 
struction, the posterior covariance Cn = {K — l)~^SnSn satisfies (12.91) . 

At this stage it suffices to know that such An and Tn exist, their finer properties play 
no role in the discussion concerning stability, but we refer the reader to HE]. 

Based on our description above, the augmented process {f7„, ..., is 

again a Markov chain. As in the previous section, we employ the notation IE„ to 
denote conditional expectation with respect to Fn- 
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3 Adaptive inflation methods 


In this section we introduce the adaptive inflation modification of EnKF, as well as the 
square root filters ETKF and EAKF. We start with the modification of EnKF, which we 
refer to EnKF-AI with the suffix standing for ‘adaptive inflation’. The nomenclature 
for the algorithms and their modifications is summarized at the end of the section, in 
Table [H 

The adaptive inflation algorithm is precisely the EnKF algorithm with the fore¬ 
cast covariance Cn replaced by an inflated forecast covariance Cn = Cn + \nl- More 
precisely, the filter ensemble {Vn^'^}k<K is governed by 



Cn 


kW = - CnH^{I + HCnH^)-\HV^^^ - 



Cn T Cn — 


— 1 ^ 

V — _ — z + 

^ ~ ’ ^n+l — w , 

k=l 

1 ^ - - 

® • 

^ k=l 


(3.1) 


The inflation parameter is chosen in such a way that it only plays a role when 
the filter is malfunctioning. In particular, when the filter is giving accurate predictions 
the inflation parameter will be zero. The inflation will be ‘triggered’ when either of two 
statistical properties of the filter, 0„ or exceeds appropriately chosen thresholds, 
indicating that filter predictions have begun to stray far from the truth. Hence we 
define 

= <^(0ru ^n) i u) = C^x{l + or y>M 2 } (3’^) 

where ip plays the role of a cut-off function, c,^ is some fixed positive constant and 
Ml, M 2 are fixed positive thresholds constant used to decide whether the filter is func¬ 
tioning properly or not. 

The first statistical quantity 0^, measures how far predicted observations are from 
actual observations: 



In a standard Kalman filter, the quantity HVn — Zn is called the innovation process 
and is used for to test accuracy of the filter. Hence the quantity On should be thought 
of as the norm of the ensemble innovation process. 

The second statistical property En measures the forecast ensemble covariance be¬ 
tween observed components and unobserved components. Using the observed-unobserved 
notation from Section 12.31 we have En = ||H„|| where Bn is the covariation between 
observed and unobserved components 


K 


Bn = 


K - 


-Y^{X^J^)-Xn)®{W-Yn). 


k=l 


Intuitively, it is important to control this cross-variation, since this avoids the situation 
whereby the filter magnifies a small error in the observed component and imposes it on 
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the unobserved component. This was a key mechanism leading to the filter instability 
found in [TO] , 

Although both statistics 0„ and can be intuitively connected to the performance 
of the filter, the real reason we choose to control these two statistics is a mathemati¬ 
cal one. If one attempts to prove boundedness results of Section 0] for the unmodified 
EnKF, it becomes clear that the ensembles cannot be stabilized unless one has adequate 
control on these two quantities. Hence it is a simple observation that an adaptive infla¬ 
tion which guarantees control of these two statistics will therein guarantee boundedness 
and stability of the corresponding adaptive hlter. 

3.1 Adaptive inflation for other ensemble based Alters 

The adaptive inflation approach can be applied to other ensemble based hlters for 
guaranteed stability with minimal changes to the original filtering method. 

Two popular modifications to EnKF are constant additive covariance inflation and 
constant multiplicative covariance inflation. With these two modifications, the analysis 
step in EnKF operates with the ensemble covariance Cn in (|2.5p replaced by 

Cn = Cn + pi , Cn = (1 + p)Cn , (3-3) 

respectively, and p > 0 here is a fixed constant. We will refer to EnKF with constant 
inflation as EnKF-Cl, with the suffix standing for ‘constant inflation’. Typically we 
will use additive rather than multiplicative inflation, but the theoretical results will 
hold for both verbatim. 

We can combine constant and adaptive inflation by simply adding both simultane¬ 
ously 

Cn = Cn -\- pi -\- Xnl 1 Cn = (1 + p)Cn + (3.4) 

respectively in the analysis step (12.5p . We will refer to this as EnKF-CAI, with the 
suffix standing for ‘constant-(-adaptive inflation’. 

In the case of the ESRFs introduced in Section 12.41 we follow the similar pattern 
of simply replacing Cn with some inflated Cn- However, this replacement only occurs 
when calculating the mean and not when calculating the posterior covariance. This is 
because Cn is the posterior ensemble covariance, its rank cannot exceed K — 1, while 
additive inflation can make the right hand side of ()2.9I) of rank d, and in most practical 
cases the model dimension d is larger than the ensemble size K. 

Hence the adaptively inflated version of ESRF is given by 

Vn = Vn- CnH^il + HCnH^)-\HVn - Zn), 

Cn = Cn- CnH^{I + HdnH^y^Hdn ■ 

where 

Cn = Cn + Xnl (ETKF-AI / EAKF-AI), 

Cn = Cn+ pi (ETKF-CI / EAKF-CI), 

Cn = Cn + pI + Xnl (ETKF-CAI / EAKF-CAI). 


(3.5) 


(3.6) 


Notice that we use the same nomenclature from Table [T] as in the case of EnKF. The 
adaptive inflation strength An defined the same as EnKF, but with 0„ now defined as 

k=l 

(k) 

since Zn no longer plays a role in the setting of ESRF. 


Acronym 

Algorithm 

Eqnation reference 

EnKF 

Ensemble Kalman filter 

(I2.5P 

EnKF-AI 

Ensemble Kalman filter w/ adaptive inflation 

eu 

EnKF-CI 

Ensemble Kalman filter w/ constant inflation 

(EJD 

EnKF-CAI 

Ensemble Kalman filter w/ constant + adaptive inflation 

(I33D 


Table 1: Acronyms for the ensemble methods. 


4 Energy principles 

In this section we show that the adaptive filtering methods introduced in Section [3] 
inherit an energy principle from the underlying model. Using the notation for the 
model introduced in Section [2l the energy principle for the model takes the following 
form: 

Assumption 4.1 (Energy principle). There exist constants j3h € (0,1), A/t > 0, sueh 
that 

E„_i|[/n|^ < (1-/3h)|Un-l|^ + A/, a.s. , 

reealling thatEn-i denotes eonditional expeetation with respect to the a-algehra J-n-i '■= 
a{Uo,Ui,... ,Un-i)- Equivalently, 

+ tr{Rh{u)) < (1 - I3h)\u\^ + Kh , 


for all for all u E 

Remark 4.2. There is a slight abuse of notation here, since we also use to 

denote expeetation eonditioned on the history of the signal-ensemble process. However 
it is elear that in the above situation the two eoineide. 

When Un is given by the discrete time formulation of an SDE as in (12.2p . it follows 
from a simple Gronwall argument that Assumption 14.11 holds as long as {'ip{u),u) < 
—/3|ttp + kh with constants I3,kh > 0. Therefore, the discrete time formulation of 
truncated Navier Stokes equation, Lorenz 63 and 96 models all satisfy Assumption 14.11 
explicit verifications can be found in Section 2.5 of [7]. 

It is natural to ask whether a filter inherits an energy principle from its underlying 
model and hence inherits a type of stability. In [7j it was shown that the non-adaptive 
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ensemble based filters (EnKF, ESRF) inherit a stronger notion of energy principle 
called the observable energy principle. In particular, if the model satisfies 

\H^h{u)\^ + tv{HRh{u)H^) < (1 - l3h)\Hu\^ + , (4.1) 

then the filter satisfies a related energy principle, which a Lyapunov function whose 
sublevel sets are only compact in the case where H is of full rank (complete observa¬ 
tions). This is clearly a much weaker result than one would hope for, since one must 
assume (14.Ih which is not thought to be true for partially observing filters and moreover 
the result is not as strong as desired since the Lyapunov function cannot be used for 
ergodicity arguments. 

We will now show that the adaptive inflation filters completely avoid these issues. In 
particular, the adaptively inflated filters inherit the usual energy principle (Assumption 
KB) and moreover, the Lyapnuov function for the ensemble has compact sub-level sets 
for any choice of observation matrix. The inflation mechanism additionally ensures that 
the observed directions of the filter HV^ never stray from the (perturbed) observations 
this is indeed the main purpose of the inflation. 

Theorem 4.3. Let denote the EnKF-AI ensemble described by (j3.ip with 

an inflation strength A„ > (p{Qn,'^n), then we have the following estimates: 

(i) The ensemble innovation satisfies the almost sure bound 

< VKma^{MuPo^c-^} 

Here po denotes the minimum eigenvalue of HqHq . 

(a) Suppose the model Un additionally satisfies an energy principle (Assumption \4-l\ ) 
and let 

£n = 4Kp^^/3fl^\\Hf\Un\^ + ^ . 

k=l 

Then there exists constant G (0,1), iA > 0 such that 

lEn-l^n < (1 — \flh)£n-l + T) . 

In particular, we have the time uniform bound sup„>oIEf’n < oo. 

Remark 4.4. It will be clear from the proof that the above result also holds verbatim 
for EnKF-CAI, introduced in (|3.4p . 

A similar energy principle also holds for the adaptive inflation modifications of 
ESRF. Note that only the second part of the result holds for the ESRF case. 

Theorem 4.5. Let denote the ETKF-AI / EAKF-AI ensemble described by 

(133]) , dMj) with any inflation strength Xn > ip{Qn,'^n), then Theorem \4 -3] (ii) holds. 

Remark 4.6. As above, it will be clear from the proof that the above result also holds 
verbatim for ETKF-CAI / EAKF-CAI introduced in (13.61) . 

We now provide the proof of Theorem 14.31 the proof Theorem 14.51 is similar and is 
therefore deferred until the appendix. 


10 


Proof for Theorem ^.5 Before proceeding, recall the ‘observed-unobserved’ notation 
from Section [2]3) 


H = {Ho,0), yW = (xW,yW), Cn = 


Cj, Bn 

Bl, Cl 


= dn = 

We start with part (i). From the assimilation equation (13.ip . by applying H to both 
sides, rearranging and applying the identity HCnH^ = HqCIHq we obtain 

_ zW = (/ + - Zll. (4.2) 

Since the minimum eigenvalue of matrix HqHq is po > 0 and HqCIHq is positive 
semi-definite, we have 

/ + HoClHl = 1 + HodlHl + ifiQn, ^n)HoHl ^ (1 + Hn))/d ■ 


Bn 

tdT 

^ r) 


Therefore, when < Mi, we have 




f T PoT(.®n^ ‘^n) 


< - Zl^\ < VkMi 


On the other hand, when > Mi, we have 








1 “1“ '^n) PO^ip^ni,^ “t“ ■^n) 


< 




PO^ip^n 

< VKpq^c-\ 


This bound and (14.2|) yield claim (i). 

We now move on to part (ii), starting by bounding the observed part. Starting from 

the expression (12.7p and applying part (i), Young’s inequality and the independence of 
(k) 

, we obtain 


E„_i|xW| 2 < 2E„_i (|XW - 

< 2po 

— 2/^0 ^K max{Mi, p^ -|- 2pQ ^||M||^E„_i|17np -|- 2pg -|- 

= 2po^iFmax{Mi,po^c“^} -h 2po ^||M||^E„_i|t/„p +4pQ^d. (4.3) 

For the unobserved part, starting from (|2.8p and applying Young’s inequality (Lemma 

EtD 


< (l+^/3,)|yi'=f+(l+2/3^-i)|B^Mo^(/+MoCfMo^)-i(MoYW-yW)|2 , (4.4) 
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where jSh G (0,1) is the one appearing in AssumDtion l4.ll For the first term on the right 
hand side, using Assumption 14.11 and the elementary inequality (1 + ^/3/i)(l — /3h) < 
(1 — we have 

E„_i(l + < (1 + < (1 - + (1 + . 


Since = ||i?n||) the second part of (14.41) is bounded by 

<KelEl\\Hof\\{I + HoC^H^)-^f 

where in the last inequality we used the bound ~ Zn'^\‘^ < KQ^. Notice that, 


as in the proof of part (i) 


/ + HoC^H^ ^ (1 + XnPo)Ig 


Hence, when 0„ > Mi or > M 2 we have 






< Po^ 


V ’ 


and when 0^ < Mi and En < M 2 , 

||(/ + MoC^Mo^)-i||0„E„ < 0„H„ < M 1 M 2 . 


In conclusion, the following holds almost surely 

\B^Hl (/ + < K\\Ho\? max{Mf M|, ^^-2} . ( 4 . 5 ) 

By combining the above estimates with (|4.4I) . the sum of (14.31) and (|4.4I) over all k 
yields the following inequality with a constant D 2 > 0 

E„_i |hW| 2 - 2A:po-i||Mof E„_i|C/„ 12 < (1 - i/3,) ^+ D 2 . (4.6) 

k=l k=l 

Then notice a multiple of Assumption 14.11 is 

AKp-^p^^\\HfW.n-i\Un? < AKp-^p^^\\Hf{l - (5h)\Un-i\^ + AKp-^p^^\\HfKh. 

The sum of the two previous two inequalities yields 

E„_iT, = + 4A:/3-VoHi - -/3/^)||M||X-l|C/„i 

k=l 

K 

= - 2i^PoiMfE„_i|C/„|2+ 4iF/3;7Vo1^HX-i|C/ni 

fc=i 

< (1 - i/3,) ^ + 4i^/3-Vo'll^li(l - mUn-i? + D 

k=l 

1 ^ 1 1 

< (1 - 2 ^") E 'll^li(l - + D = {1- -MSn, 

k=l 

with D := Pq^ ||i/|piF, + D 2 and using the fact that (1 — /?,) < (1 — \l3h)^ ■ D 


+ D 2 , 
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5 Geometric ergodicity 

A stochastic process is geometrically ergodic if it converges to a unique invariant mea¬ 
sure geometrically fast. In this section we will show that the signal-ensemble process 
(for each filter introduced in this article) is geometrically ergodic. In particular, if P 
denotes the Markov transition kernel for the signal ensemble process, then we will show 
that there exists a constant 7 G (0,1) such that 

\\P^^i-P^u\\TV<C^,ur , (5.1) 

where /x, v are two arbitrary initial probability distributions, is a time uniform 
constant that depends on /i,zx, and || • \\tv denotes the total variation norm. 

Geometric ergodicity is a notion of stability for the filter. In particular it implies 
that discrepancies in the initialization of the filter will dissipate exponentially quickly. 
For the linear Kalman filter and optimal filters, it is known that under mild assumptions 
geometric ergodicity is extended from model to filter [221 [23] . 

In theorems 5.5, 5.6 and 5.8 of [7], the authors proved that if the system noise Cn is 
non-degenerate and the signal-ensemble process has an energy principle, then EnKF, 
ETKF and a natural version of EAKF will be geometrically ergodic. In particular, non 
degeneracy for the system noise Cn means the following. 

Assumption 5.1 (Nondegenerate system noise). For any constants > 0, there 

is a constant a > 0 such that 

P(Cn G • \Un-i =u)> a\R^{-) 

for all |tt| < Ri, where \R^{dx) is the Lebesgue measure ofW^ restricted to {u : |u| < 
R 2 ]. 

Assumption 15.11 holds for many practical examples. When Un is produced by time 
discretization of an SDE (12.21) . it suffices to require S being nonsingular, see Appendix 
B in [7] for a detailed discussion. 

In contrast with the ergodicity results of [7j , we no longer require the assumption of 
a Lyapunov function with compact sub-level sets, since this is guaranteed by Theorems 
K2\ and 1131 This is a vast improvement as in the non-adaptive case such Lyapunov 
functions are only known to exist in the (essentially) fully observed scenario. 

Before stating the result, recall that Mi, M 2 > 0 are the thresholds used to decide 
when to trigger the adaptive inflation mechanism. For EnKF-AI and ETKF-AI, the 
thresholds have no constraints besides being positive, but for EAKF-AI we require that 
the constants are sufficiently large. 

Theorem 5.2. If the system noise in the signal process is non-degenerate (Assump- 
tion \5.1\) then the signal-ensemble processes {Un, Vn^\ ..., Vn^^) generated by EnKF-AI 
(with any Mi, M 2 > 0), ETKF-AI (with any Mi, M 2 > 0) or a version of EAKF-AI 
(with Ml, M 2 sufficiently large) are geometrically ergodic. 

Remark 5.3. The same result holds verbatim for adaptive filters with additional con¬ 
stant covariance inflation (EnKF-CAI,ETKF-CAI,EAKF-CAI), which are defined in 
Section I3.il 

As a matter of fact, the above result is an immediate corollary of [Tj- The reasoning 
will become clear after we review the major steps in [7|. 
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5.1 A controllability framework 

The central piece of the arguments is a result of Markov chain theory [23]. Here we 
use a simple adaptation of the form given in |25t Theorem 2.3]. 

Theorem 5.4. Let X„ he a Markov chain in a space E such that 

1. (Lyapunov function) There is a function £ : E R"*" with compact sub-level sets 
and such that 

E„_iT(X„) < (l-/3)£:(X„_i)+X, 
for some (3 G (0,1) and K > 0 . 

2. (Minorization condition) Let Cr = {x : £{x) < R}. For all R > 0, there exists 
a probability measure n on E with v{Cr) = 1 and a constant r] > 0 such that for 
any measurable set A C E 

P(X„ G A|X„_i = x) > r]u{A) 

for all x G Cr. 

Then there is a unique invariant measure vr and constants r G (0,1), k > 0 such that 


'(Xn G •) — -kWtv < ( 1 + / £{x)pi{dx) 


Theorems 14.31 and 14.51 have already provided the Lyapunov function. All that 
remains is to show the minorization condition. 

To obtain minorization, we will exploit the update structure of the ensemble filters, 
precisely as was done in [7j. Notice that for all ensemble filters, the signal-ensemble 
process X„ := ([/„, Vn^\ ..., Vn^'^) is a Markov chain taking values in A := R'^ x 
The evolution of X„ is described by the composition of two maps. The first is a random 
map from A to a signal-forecast-observation space T, described by a Markov kernel 
<f> : A X B{y) [0,1]. The second is a deterministic map T : T ^ A, which combines 
the forecast with the observed data to produce the updated posterior ensemble. The 
details of these maps, as well as the definition of the intermediate space T, differs 
between EnKF, ETKF, EAKF and their adaptive counterparts. 

For example, in EnKF, the intermediate space is T := R'^ x R'^^^ x R*?^^ and the 
random mapping is 


{Un-I, ..., V^^l) ^ := {Un, .., fW, • ■ ■, 


riK) 


The deterministic map T is given by the Kalman update 


TiUn, w ,..., .., zf'>) = {Un, r(i),...,rW), 


where 


p(fc) ^ y{k) _ 

K 


c = 


1 


K - 


'^{yik) y) 


(5.2) 


k=l 


For EnKF-AI, we have the same formulas, but with the exception 
K 


C = 


1 


K - 


_ - F) ® (fW - F) + c^0„(l + S„)le„>Mi or 


/ . 


k=l 
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and 


0 


n 



1 

K - 1 


K _ 


k=l 



The corresponding formulas for ETKF, EAKF (and their adaptive counterparts) can 
be similarly derived and for concreteness are given in Appendix O Notice that for 
values of close to the origin, the adaptive inflation term will vanish and hence the 
deterministic map E for EnKF-AI coincides with that of EnKF near the origin. Similar 
statements hold for ETKF-AI and EAKF-AI. 

Given this formulation, it suffices to show that the push-forward kernel r*$(x, •) = 
‘h(x,r~^(-)) satisfies the minorization condition. It is easy to see that, given the 
assumptions on the noise, the kernel d>(x, •) has a density with respect to Lebesgue 
measure, so we simply need to show that the pushforward inherits the density property 
from <I>. This can be understood as the controllability of the map T. To achieve this, 
we use the following simple fact, which is lemma 5.4 in [7]. 

Lemma 5.5. Let ^ he a Markov transition kernel from ML xR™ with a Lebesgue 

density p{x,y) = p{x, {yi,y 2 )) and let T : R” x R™ ^ R”. Given a compact set C, 
suppose that there is a point y* = {yl,y 2 ) £ x and /3 > 0 such that 

1. Reachable from all x £ C, that is the density function p{x,y) > /3 for y around 

y*, 

2. T is controllable around y*, that is T is near y* and det(Py^r)|j^* > 0 . 

Then there is a 6 > 0 and a neighborhood Oi ofT{y*) such that for all x G C 


r$(x,-) >5Aoi(-) 


where Aoi is the Lebesgue measure restricted to the set Oi. In other words, the mi¬ 
norization condition holds for the transition kernel r*<k. 


5.2 Application to ensemble filters with adaptive infiation 

To apply Lemma 15.51 to the proof of Theorem 15.21 we will use the variables 

X = ([/„_!, 

yi = ([/„, E«,...,eW) (5.3) 

2/2 = (Z„,zW,...,zW) 

The choice of the intermediate point y* = {y\,y 2 ) can be quite delicate and can sim¬ 
plify the non-degeneracy condition considerably. In particular, we prove the following 
proposition in [7]. 

Proposition 5.6. Let T he the Markov kernel and deterministic map used to define 
EnKF/ETKF/EAKF, given by the formulas in AppendixUl. If the system noise is non¬ 
degenerate (Assumvtion \5. 1\) then 4>,r satisfies both conditions o f Lemma 1 5 . 51 with the 
following choice of intermediate point y*: 

1. EnKF: y* has all its components being at the origin. 
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2. ETKF: y* has all its components being at the origin. 

3. EAKF: has all its components being at the origin, U* is at the origin, while 

..., is the d x K matrix given by 


1 



if d<K -1 


1 Ho 
0 0 


if d>K -1 


where 1 is the dxl vector of 1 ’s and Hq is the rxr matrix with r = minji^r — 1, d}: 


1 1 ... 1 -rl 


1 -2 0 
-10 0 


0 

0 

0 

0 


Note in this case, a version of EAKF is picked so T is near y*. 


With this propsition at hand, Theorem 15.21 becomes an immediate corollary. 


Proof of Theorem \5.2[ For each of EnKF-AI, ETKF-AI, EAKF-AI verification of the 
reachability condition in Lemma [53] is identical to that of EnKF, ETKF, EAKF which 
is trivial and is done in [3 Proof of theorems 5.5, 5.6, 5.8] respectively. Hence it 
suffices to check the non-degeneracy of F in each of the three cases. For EnKF-AI and 
ETKF-AI we pick y* = ( 2 / 1 , 2 / 2 ) to be the origin. Hence for thresholds Mi, M 2 > 0, 
there exist a neighborhood of the origin such that the map F coincides with that of 
EnKF, ETKF respectively, since the inflation mechanism is not triggered within this 
neighborhood. Hence the non-degeneracy follows immediately from Proposition 15.61 
For EAKF-AI, as detailed above, the intermediate point y* is not chosen at the origin 
but rather at the point specified in Proposition 15.61 Nevertheless, since this point is 
fixed, one can pick thresholds Mi, M 2 sufficiently large so that for some neighborhod 
of y* the inflation mechanism is not triggered and hence F agrees with that of EAKF. 
Hence the non-degeneracy follows from Proposition 15.61 □ 


6 Thresholds from elementary benchmarks 

Finally we discuss the choice of the cutoff thresholds Mi and M 2 in the formulation 
(|32D of (p. While Theorems 14.31 guarantees filter stability for any Mi, M 2 > 0, certain 
values produce better filter performance. Intuitively, these thresholds should be chosen 
to differentiate malfunctioning forecast ensembles from properly working ones. 

6.1 An elementary benchmark 

Here we devise an elementary benchmark of accuracy that should be surpassed by any 
properly functioning filter. At each observation step, the benchmark estimate is given 
by the stationary state Un conditioned on the observation This can be computed 
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using Bayes’ formula and only requires access to the invariant measure vr of the model. 
For instance we can compute the conditional expectation via 


E(C/„|Z„) = 


f X7r(x) exp(—— IIxl^)dx 


( 6 . 1 ) 


f 7r(x) exp(— 1 |Z„ — ffxl^)dx 

The benchmark for accuracy is given by the mean-square error of the above estimator; 

Errors := E\Un - E{Un\Zn)\^ = E\Un\‘^ - E\E{Un\Zn)\‘^ . 

A properly functioning forecast ensemble should be expected to perform better than 
this estimator. For instance, one should expect a properly functioning filter to satisfy 


1 

K 


K 




k=l 


< Error. 


Consequently, one should expect that 
K 

K 


1 ^ 

Eel = ^ - Un) +Cn + < 4 := Errors + 2d . 


k=l 


Moreover using the inequalities ||a® 6 || < |o||5|, |x||y| < + y^) and Lemma fA. 2 1 we 

expect 


K 


E 






< 


1 


K 


k=l 
K 


2K -2 




k=l 


< -^E|Ei*’) - Unl"^ < Ms := ^ Errors. 


2K-2 


k=l 


2K-2 


In summary, we expect that a properly functioning filter should satisfy 

E0^ < 0-0 := ||iL|pErrorA + 2d, ES„ < Ms ;= Errors . ( 6 . 2 ) 

Since the integrals appearing in the Bayesian formula (j 6 .ip are often difficult to com¬ 
pute, it is more useful to have an estimator that only relies on low order statistics of 
vr. For this purpose it is natural to introduce a Gaussian approximation of tt, replacing 
vr with TT = M(EUn, cov(Un)) where the mean and covariance are the true statistics 
computed from vr. The condtional distribution of vf given the observation is now a 
Gaussian measure and can be computed exactly. In particular we have that 

E{Un\Zn) = EUn - COv{Un)H^{I + COv{Un)H)-\HEUn - Z^). 

The covariance of the error = Un — E{Un\Zn) is given by the Kalman prior-posterior 

covariance relation 

cov(ri^) = COv{Un) — COv{Un)H'^ {I + H COY COY {Un), (6.3) 

which can be explicit computed. The mean square error is given by 

Errors = E|ri^p = tr(cov(r/f)). 
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In the case of a linear model, this Gaussian approximation is of course exact. 

Given the above discussion, an obvious option for the thresholds is 

Ml = cT©, M 2 = Ms. (6.4) 

we call this aggressive thresholding. In this situation the adaptive inflation is triggered 
as soon as the hlter statistics 0^, exceed the mean values computed from the crude 
estimators. A less aggressive strategy is to obtain Mi, M 2 by hypothesis testing, which 
we do not pursue here. 

Remark 6.1. To avoid wasting resources on computing these benchmarks, one could 
instead use the trivial estimator Un = 0. In this case, one can compute 

nUn-Q\'^ =nUn? <l3^^Kh. 

A properly functioning EnKF should have a smaller mean square error in its forecast. 
And as a consequence, we can replace term Error a in the bounds above with 
which leads to 

EQl < + 2d, EH, < 

7 Numerical results 


In this section, we will numerically validate our results through the following five mode 
Lorenz 96 model: 


d 

— Xi = Xi- 


i{xi+i - Xi- 2 ) -Xi + F, 


i = 1,...,5. 


(7.1) 


Here the indices are interpreted in modulo 5 sense. The truth will be given by Un = 
x{nh) for some observation time step h. Previous studies from |26l |9l 0] have shown 
that different choice of forcing F produces very different strength of turbulence. Here 
we will pick forcing T = 4,8 and 16, and see that these choices create very different 
filtering results. 

Throughout the section we assume that only one component xi is observed, hence 
we have the observation matrix H = [1,0, 0,0, 0]. The observational noise is chosen to 
be distributed as AA(0,0.01). We apply a selection of ensemble assimilation methods 
to this signal-observation system with 6 ensemble members, where the initial positions 
are drawn from a Gaussian approximation of the equilibrium measure. For the sake of 
simplicity, we integrate (17.ip using an explicit Euler scheme with time step A = 10“^. 
For most ensemble based filtering methods, the explicit Euler scheme is bad choice of 
integrator for (17.ip . since the stiffness of the equation will lead to numerical blow-up 
for sufficiently large initial conditions. As we shall see, this leads to prevalence of 
catastrophic filter divergence in most ensemble based filtering methods. The adaptive 
inflation method on the other hand avoids this issue and exhibits high filtering skill. 

Here we compare the performance of four ensemble based filters: EnKF, EnKF-AI, 
EnKF-GI and EnKF-CAI, which have been defined in Sections [2] and [3l The constant 
inflation is always taken to be additive. 

Following the methodology introduced in Section [6l we construct the equilibrium 
measure of the Lorenz 96 system using data from a long time (T = 10^) simulation, 
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which leads to the performance benchmarks computed by (16.311 and (16.21) . The thresh¬ 
olds Ml and M 2 are obtained through the aggressive strategy (j6.4l) . using the Gaussian 

approximation of the stationary measure tt. In order for the comparison to be unbi- 

(k) 

ased, we use the same realizations of observation noise and perturbations ' in all 
four filters. The statistics are collected for N = 100 independent trials, where each 
runs a total time length of T = 100. To eliminate transients, we will only record data 
from the period 50 < T < 100. 

The performances of the four filters are measured and compared from three per¬ 
spectives: 

• The most important question here is how often does catastrophic filter diver¬ 
gence appear in standard (non-adaptive) filters, like EnKF and EnKF-CI, and 
does the adaptive inflation in EnKF-AI and EnKF-CAI prevent it. Catastrophic 
filter divergence can be identified by checking whether the ensemble at the final 
time T = 100 takes value “NaN”, representing machine infinity, in any of its 
components. 


• The overall accuracy of the filters is compared using the average of the root mean 
square error (RMSE) over N = 100 trials among the second half of the time 
interval, where 

2 - 

RMSE= - ^ \Vn-Un\^. 

n=Tl2h 


Another statistics that is useful in judging the accuracy is the averaged pattern 
correlation: 


Cor 


2 ^ (Vn U) 


where U denotes the climatological average of the true model (the mean of the 
invariant measure) and Vn as usual denotes the average of the filter ensemble. 
We will also plot the of the posterior error, \Vn — Un\, for one realization of a 
typical trajectory for the four filters. 

The following experiments have also been carried out with both ESRF methods. 
All these methods have performance very similar to EnKF based method, as well as 
the effect of adaptive inflation over them. We do not present these results for brevity. 


7.1 Accuracy for different turbulent regimes 

In this section we compare the performance of the four filters across the three turbulent 
regimes. In each experiment we use observation interval h = 0.05 and (when required) 
constant additive inflation strength p = 0.1. 

7.1.1 Weak turbulence 

Here we compare the performance of the filters in a weakly turbulent regime, by setting 
F = A. The performances of the four filters are presented in Table [2l where we compare 
frequency of catastrophic filter divergence, RMSE and pattern correlation. In Figure 
EH we compare the posterior error for the four filters for one typical trajectory. 
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We found that, among 100 trials, none of the filters experience catastrophic filter 
divergence. The adaptive inflation in EnKF-AI has been triggered in 30 trials, and on 
average 1.96 times in each trial. In EnKF-CAI, the adaptive inflation has been turned 
on in 9 trials, and only once in each trial. Hence the typical realizations of EnKF and 
EnKF-AI are identical and similarly for EnKF-CI and EnKF-CAI, since the adaptive 
inflation is seldom triggered. 

In terms of accuracy, all four filters perform quite well. For each filter the RMSE is 
significantly lower than the benchmark RMSE, which is 3.64 and the pattern correlation 
is very high. Unsurprisingly, the filters with constant covariance inflation (EnKF-CI 
and EnKF-CAI) have more skill than EnKF, with lower RMSE and higher pattern 
correlation. Somewhat more surprising is that EnKF-AI also performs better than 
EnKF, both in RMSE and pattern correlation, even though the adaptive covariance 
inflation has only been switched on in 30 of the 100 trials. This indicates that the 
performance of EnKF must be quite bad in those 30 and the adaptive inflation is 
needed to ensure the accuracy. Of course, in this regime one could also achieve this 
with non-adaptive constant inflation. 


Filter 

EnKF 

EnKF-AI 

EnKF-CI 

EnKF-CAI 

Cata. Div. 

0% 

0% 

0% 

0% 

RMSE 

0.89 

0.54 

0.22 

0.22 

Pattern Cor. 

0.91 

0.96 

0.98 

0.98 


Table 2: In the weak turbulence regime {F = 4), for the four algorithms we compare 
frequency of catastrophic Alter divergence (Cata. Div.), RMSE, and pattern correlation 
(Pattern Cor.) The benchmark RMSE is 3.25. 


We now briefly describe how the thresholds and benchmarks were calculated. The 
equilibrium mean on each mode is Xi ~ 1.22 and the variance is var(xi) ~ 3.38. 
Using Kalman one time assimilation, the posterior variance of xi will be reduced to 
0.01, while the benchmark RMSE will be 3.25. From these statistics, we compute 
Ml = 0-0 = 32.5, Ms = Ms = 6.2. 

7.1.2 Moderate turbulence 

Here we compare the performance of the Alters in a moderately turbulent regime, by 
setting F = 8. The performances of the four filters are presented in Tabled where we 
compare frequency of catastrophic Alter divergence, RMSE and pattern correlation. In 
Figure [721 we compare the posterior error for the four filters for one typical trajectory. 

Among 100 trials, catastrophic filter divergence takes place 12 times for EnKF, but 
never takes place for the other three. The adaptive inflation in EnKF-AI has been 
turned on in 96 trials, while on average 17.5 times in each of these trials. As for EnKF- 
CAI, the adaptive inflation has been turned on in 20 trials, while on average 6.05 times 
each of these trials. 

In terms of accuracy, EnKF will not be analyzed due to the frequent occurrence of 
catastrophic Alter divergence. The EnKF-AI avoids the issue of catastrophic filter di¬ 
vergence, however it has very little skill, since the RMSE is higher than the benchmark 
RMSE 7.02 and the pattern correlation is quite low. However, the filters with constant 
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Figure 7.1: We compare the posterior error for one realization of the four filters in the weakly 
turbulent (F = 4) regime. Since the adaptive inflation is not triggered, trajectories collapse 
into two groups (those with constant inflation and those without). The benchmark RMSE 
is included (black dash) to emphasize the accuracy of the hlters. 


inflation (EnKF-CI and EnKF-CAI), behave almost identically, display significantly 
better skill, beating the benchmark RMSE and having high pattern correlation. 


Filter 

EnKF 

EnKF-AI 

EnKF-CI 

EnKF-CAI 

Cata. Div. 

12% 

0% 

0% 

0% 

RMSE 

NaN 

8.6 

3.61 

3.57 

Pattern Cor. 

NaN 

0.55 

0.89 

0.89 


Table 3: In the moderate turbulence regime (F = 8), for the four algorithms we compare 
frequency of catastrophic hlter divergence (Cata. Div.), RMSE, and pattern correlation 
(Pattern Cor.) The benchmark RMSE is 7.02. 


We now briefly describe how the thresholds and benchmarks were calculated. The 
equilibrium mean on each mode is Xi ~ 2.28 and the variance is var(xj) ~ 12.6. 
Using Kalman one time assimilation, the posterior variance of xi will be reduced to 
0.01, while the benchmark RMSE will be 7.02. From these statistics, we compute 
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Figure 7.2: We compare the posterior error for one realization of the four hlters in the 
moderately turbulent {F = 8) regime. The EnKF and EnKF-AI are not performing well, as 
they are frequently exceeding the benchmark (black dashed). Adding constant covariance 
inflation signihcantly improves the performance. 


Ml = (Te = 69.56 and M 2 = M= = 28.8. 


7.1.3 Strong turbulence 

Here we compare the performance of the filters in a highly turbulent regime, by setting 
F = 16. The performances of the four filters are presented in Table 01 where we 
compare frequency of catastrophic filter divergence, RMSE and pattern correlation. In 
Figure ITAI we compare the posterior error for the four filters for one typical trajectory. 

In this highly turbulent regime, both non-adaptive filters, EnKF and EnKF-CI 
display frequent catastrophic filter divergence, making them practically useless. Among 
100 trials, catastrophic filter divergence occurs in every trial of EnKF and in 18 trials of 
EnKF-CI. Both EnKF-AI and EnKF-CAI experience no catastrophic filter divergence. 
The adaptive inflation has been turned on in all trials of EnKF-AI, occurring 97.53 
times on average per trial. The adaptive inflation has been turned on for EnKF-CAI 
in 80 trials, occurring 18.7 times on average per trial. 

In terms of accuracy, the EnKF-AI exhibits very little skill, with RMSE significantly 
higher than the benchmark RMSE 12.93 and with low pattern correlation. The EnKF- 


22 








































































CAI however does slightly beat the benchmark RMSE and displays moderate pattern 
correlation. In Section 17.21 below, we show that the performance of EnKF-CAI is 
sharply improved by a simple tuning of parameters. 


Filter 

EnKF 

EnKF-AI 

EnKF-CI 

EnKF-CAI 

Cata. Div. 

100% 

0% 

18% 

0% 

RMSE 

NaN 

24.48 

NaN 

11.91 

Pattern Cor. 

NaN 

0.23 

NaN 

0.69 


Table 4: In the highly turbulence regime {F = 16), for the four algorithms we compare 
frequency of catastrophic hlter divergence (Cata. Div.), RMSE, and pattern correlation 
(Pattern Cor.) The benchmark RMSE is 12.93. 




Figure 7.3: We compare the posterior error for one realization of the four hlters in the 
strongly turbulent {F = 16) regime. The hlter is very unstable, EnKF (red, upper panel) 
explodes to machine inhnity, and EnKF-CI (red, lower panel) does the same with signihcant 
probability. The adaptive inhation mechanism, when applied to these two hlters, is triggered 
frequently, preventing the ensemble from exploding. The performance of EnKF-AI (blue, 
upper panel) is worse than the benchmark (black dash). The performance of EnKF-CAI 
(blue, lower panel), is much better and outperforms the benchmark. 
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We now briefly describe how the thresholds and benchmarks were calculated. The 
equilibrium mean on each mode is Xi ~ 3.1 and the variance is var(xi) ~ 40.6. Using 
Kalman one time assimilation, the posterior variance of xi will be reduced to 0.01, 
while the benchmark RMSE will be 12.93. From these statistics, we compute Mi = 
iT0 = 127.6 and M2 = Ms = 81.4. 
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7.2 Dependence on constant inflation strength and obser¬ 
vation time 


Based on previous discussions, the performance of EnKF-CI and EnKF-CAI are sat¬ 
isfactory in the weakly and moderately turbulent regimes, but the strongly turbulent 
regime poses a significant challenge. In this section, we study how constant covariance 
inflation strength p and observation time h affect the performance of EnKF-CI and 
EnKF-CAI in the E = 16 regime. 

In Table [5] we compare the performance of EnKF-CAI for different choices of con¬ 
stant inflation strength /?, including RMSE and pattern correlation. We also include 
the frequency of catastrophic filter divergence in EnKF-CI to emphasize the impact of 
adaptive inflation. 

As constant inflation strength is decreased, the frequency of catastrophic hlter 
divergence increases. This agrees with the notion that removing inflation simply results 
in EnKF, which in Table [5] was shown to experience catastophic filter divergence in 
100 % of trials. 

In terms of accuracy, the EnKF-CAI can be optimized by choosing a moderately 
small constant inflation strength. In particular, when p is in the 0.01 — 0.02 range 
the RMSE is significantly smaller than the benchmark RMSE 12.93 and the pattern 
correlation is relatively high. Surprisingly, the optimal choice of p for EnKF-CAI is 
very sub-optimal for EnKF-CI, with around half of all trials experiencing catastrophic 
filter divergence. This emphasizes the importance of creating the right balance between 
adaptive inflation and constant inflation. 


p 

1 

0.5 

0.2 

0.1 

0.05 

0.02 

0.01 

0.005 

EnKF-CI Cata. Div. 

3% 

8% 

18% 

18% 

28% 

42% 

57% 

75% 

EnKF-CAI RMSE 

13.05 

13.62 

13.43 

11.91 

8.82 

8.51 

9.3 

10.51 

EnKF-CAI Cor. 

0.64 

0.65 

0.66 

0.69 

0.70 

0.70 

0.75 

0.70 


Table 5: In the F = 16 regime we display for different choice of constant inflation strength 
p, the frequency of catastrophic hlter divergence for EnKF-CI (EnKF-CI Cata. Div.), the 
RMSE for EnKF-CAI and the pattern correlation for EnKF-CAI (EnKF-CAI Cor). Each 
statistic is generated from 100 independent trials. The benchmark RMSE is 12.93. 

In Table [ 6 ] we display the results of a similar experiment, now varying the observa¬ 
tion interval h keeping a fixed inflation strength p = 0.1. Again we display frequency 
of catastrophic filter divergence for EnKF-CI and the and pattern correlation RMSE. 

Each statistic is generated from 100 independent trials. 

We found that the frequnecy of catastophic filter divergence is not monotonic in 
h, but rather displays a peak at h = 0.1. This agrees with the hndings of [9]. We 
also hnd that the accuracy of EnKF-CAI is not monotonic in h but rather, there is an 
optimal choice around h = 0.1. As with the inflation strength experiment, the optimal 
parameter choice for the accuracy of EnKF-CAI occurs when EnKF-CI is exhibiting 
frequent catastrophic filter divergence. The fact that filter performance decreased with 
decreasing h can be attributed to a violation of filter observability conditions when the 
time-step becomes too small [ 1 ]. 
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h 

0.01 

0.02 

0.05 

0.1 

0.2 

0.5 

EnKF-CI Cata. Div. 

0% 

1% 

18% 

25% 

5% 

0% 

EnKF-CAI RMSE 

25.75 

20.71 

11.91 

6.43 

14.09 

14.80 

EnKF-CAI Cor. 

0.31 

0.37 

0.69 

0.64 

0.50 

0.36 


Table 6: In the F = 16 regime we display for different choice of observation interval h, the 
frequency of catastrophic hlter divergence for EnKF-CI (EnKF-CI Cata. Div.), the RMSE 
for EnKF-CAI and the pattern correlation for EnKF-CAI (EnKF-CAI Cor). Each statistic 
is generated from 100 independent trials. The benchmark RMSE is 12.93. 


7.3 Dependence on the numerical integrator 

In this section, we address the question of whether the catastrophic filter divergence can 
be avoided by simply using a more stable integrator in EnKF and EnKF-CI. We will see 
that with implicit methods and variable time step methods, one can avoid catastrophic 
filter divergence, but at prohibitive computational cost. Moreover, although they avoid 
catastrophic filter divergence, they do not avoid ordinary filter divergence; in particular 
we will see that EnKF and EnKF-CI with implicit / variable time-steps do not exhibit 
any filtering skill, in contrast with EnKF-AI and EnKF-CAI. This corroborates the 
findings of [8l Chapter 2], where it was shown that numerical model error from implicit 
methods leads to reduced filtering skill. 

In Table [TJ we display the results from 100 trials, comparing the performance of 
EnKF-AI with explicit Euler (step size A = 10“^) against EnKF with (respectively) 
explicit Euler (A = 10““^), standard 4-th order Runge-Kutta (A = 2.5 x 10“^), ode45 
and implicit Euler (A = 10“^). ode45 is a native MATLAB integrator, which is based 
on a Runge-Kutta (4, 5) method with variable step size. The implicit equation in 
implicit Euler is solved using the native MATLAB function f slove. The average time 
is also shown to indicate the computational expense associated with each method. 

We see that ode45 and implicit Euler are the only methods that do not exhibit 
catastrophic filter divergence, but each trial is enormously expensive when compared 
to the explicit Euler based methods. Moreover, ode45 and Implicit Euler display zero 
filtering skill, as judged by the RMSE and pattern correlation (the benchmark RMSE 
is 12.93). As discussed earlier, EnKF-AI also exhibits little skill in this regime, but 
we shall see that the skill is dramatically improved by introducing additional constant 
inflation. Unsurprisingly, the explicit Euler based methods are much faster than both 
the higher order and the implicit methods and the EnKF-AI method is only fractionally 
more expensive than EnKF. 

In Table [8] we display the results from a similar experiment, now comparing EnKF- 
CAI (with explicit Euler) against EnKF-CI with the same four choices of integrator. 
Catastrophic filter divergence is overall less prevalent, with ode45 and implicit Euler 
once again exhibiting no catastrophic filter divergence. In contrast with the previous 
experiment, EnKF-CAI exhibits good filtering skill as judged by the RMSE and pattern 
correlation (once again, the benchmark RMSE is 12.93). EnKF-CI with ode45 exhibits 
less filtering skill than EnKF-CAI, with RMSE not beating the benchmark and EnKF- 
CI with implicit Euler exhibits even less skill. As in the first experiment, the adaptive 
inflation is more accurate and far cheaper than the higher order and implicit methods. 
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EnKF-AI 

EnKF 

Integrator 

Explicit 

Explicit 

RK4 

ode45 

Implicit 

Cata. Div. 

0% 

100% 

92% 

0% 

0% 

RMSE 

24.48 

NaN 

NaN 

44.66 

20.38 

Pattern Cor. 

0.23 

NaN 

NaN 

0.22 

0.11 

Avg. Time 

2.49 

2.39 

19.53 

41.00 

662.48 


Table 7: The performance and time cost of EnKF-AI (explicit Euler, A = 10“'^)) against 
EnKF in the F = 16 regime. EnKF is implemented by explicit Euler (A = 10““^), RK4 
(A = 2.5 X 10“^), ode45 and implicit Euler (A = 10“^). The benchmark RMSE is 12.93. 
The time cost is measured in seconds. 



EnKF-CAI 

EnKF-CI 

Integrator 

Explicit 

Explicit 

RK4 

ode45 

Implicit 

Cata. Div. 

0% 

18% 

5% 

0% 

0% 

RMSE 

11.91 

NaN 

NaN 

14.95 

16.29 

Pattern Cor. 

0.69 

NaN 

NaN 

0.71 

0.42 

Avg. Time 

2.65 

2.34 

21.48 

37.02 

650.31 


Table 8: The performance and time cost of EnKF-CAI (explicit Euler, A = 10“^)) against 
EnKF-CI in the F = 16 regime. EnKF is implemented by explicit Euler (A = 10““^), RK4 
(A = 2.5 X 10“^), ode45 and implicit Euler (A = 10“^). The benchmark RMSE is 12.93. 
The time cost is measured in seconds. 

7.4 The distribution of and 

The filter ensemble statistics 0,^ and defined in Section [3] are used to determine 
when the adaptive inflation of EnKF-AI should be triggered. In particular, the mech¬ 
anism is triggered whenever either of the two statistics passes an appropriately chosen 
threshold. In this section we look at the distribution of these statistics in order to 
evaluate the choice of thresholds introduced in Section [H We have collected the data 
of On and from the EnKF-AI scheme from all the analysis steps. In Table [9] we 
compare the distributions of 0^ and in the three different turbulent regimes and 
compare with the choice of threshold. 

It is clear from the table that stronger turbulence significantly increase the size of 
both 0n and The thresholds Mi and M 2 are always larger than the average of 0,^ 
and but the ratios Mi/(0„) and M 2 /(H„) significantly decreases as the turbulence 
gets stronger. As a result, the percentage of 0,^ and that are beyond the thresholds 
significantly increases. Connecting this observation with the performance of EnKF-AI, 
which is good in F = 4, slightly off the benchmark in F = 8, and very bad in F = 16, 
we can conclude that whether 0„ and are beyond the thresholds is a good indicator 
of the malfunctioning. As intended, the adaptive inflation is rarely triggered when the 
filter operates properly, but frequently triggered when the situation is chaotic. 

The histograms of 0n and in three regimes are presented in Figure 17.41 0„ 
has a Gaussian like tail in the weak turbulence regime and exponential like in strong 
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F 

4 

8 

16 

Ml 

32.5 

69.6 

127.6 

Average 0„ 

2.95 

20.28 

70 

P(0n > Ml) 

0.0285% 

3.025% 

9.73% 

M 2 

6.2 

28.8 

81.4 

Average 

0.001 

0.0437 

0.22 

P(S. > M 2 ) 

0% 

0% 

0% 


Table 9: For the three turbulent regimes {F = 4,8,16), we list the thresholds Mi and M 2 
chosen according to the method of Section O We also list the average of 0„ and over 
100 independent trials and also give the probability that the statistics pass their respective 
threshold. 


turbulence. The tail of is much heavier and polynomial-like. The thresholds we 
compute from the benchmarks are fairly large comparing with the distribution. These 
two posterior observations suggest that using even an aggressive threshold strategy 
as 16.41 the outliers are relatively scarce so adaptive inflation does not hamper the 
performance of the filters. 

It is worth noting that is usually a small number and very rarely exceeds the 
threshold M 2 . This suggests that the filter performance would be unaffected if the 
adaptive inflation was changed to so that only 0„ is used to trigger the 

inflation. We have tested this simpler adaptive inflation mechanism with both EnKF- 
AI and EnKF-CAI in the setting of the previous subsections, and the performances 
are much as the same. From a theoretical perspective however, both statistics must be 
included to obtain a rigorously stable filter. 

7.5 Discussion of Experiments 

In Section 17.11 the numerical simulations convincingly demonstrate that catastrophic 
filter divergence is prevalent in EnKF and EnKF-CI, particularly when using an un¬ 
stable integrator in the forecast step. By comparing with EnKF-AI and EnKF-CAI, 
adaptive inflation is shown to successfully avoid catastrophic filter divergence, without 
sacrificing any filtering skill. In the weakly turbulent regime, the adaptive inflation 
mechanism is rarely turned on and EnKF-AI performs similarly to EnKF and like¬ 
wise for the constant inflation counterparts. In more turbulent regimes, the adaptive 
inflation is regularly triggered and the adaptive inflation filters display filtering skill, 
typically beating the RMSE benchmark. In Section 17.21 it is also shown that this skill 
can be optimized by tuning constant inflation strength. 

We saw in Section ES that catastrophic filter divergence can be avoided by using 
implicit or variable time step integrators, but at the expense of prohibitive computa¬ 
tional cost and the loss of filtering skill. Indeed, one would be better off using the naive 
estimator developed for the RMSE benchmark. 

The adaptive inflation method permits for the effective use of very cheap integrators 
such as explicit Euler in the forecast step. As evidenced by EnKF, it is impossible to use 
such integrators in turbulent models, due to prevalence of catastrophic filter divergence. 
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Figure 7.4: Histograms of 0„ (upper panel) and (lower panel) in F = 4, 8,16 regimes. 
The thresholds Mi are marked as black dashes. The thresholds M 2 lie out side the x-axis 
for Note the logarithmic scale on the vertical axis. 


The adaptive inflation completely avoids this issue. 

In Section 17.41 we provide the distribution analysis of 0^ and which indicates 
the thresholds are achieving our purpose. We also saw that the threshold for was 
never triggered, indicating that a simpler adaptive inflation only using the statistic 
0n will have practically the same performance. This has been validated, but is not 
presented here. 


8 Conclusion and Discussion 

In this article we have proposed a novel modification of ensemble based methods using 
an adaptive covariance inflation. This modified method resolves common issues with 
ensemble based methods, which were raised in the introduction, namely that there is 
a lack of understanding of the stability properties of EnKF, with catastrophic filter 
divergence drawing stark attention to this fact. The EnKF-AI and related filtering 
methods (ETKF-AI, EAKF-AI, EnKF-CAI, etc.) resolve this stability issue entirely. 

In particular we can develop a complete stability theory for these filters. In Section 
0] we have proved time-uniform mean-square upper bounds for EnKF-AI and related 
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filters, using a Lyapunov argument. This is a considerable improvement on the known 
results for EnKF [7] since the observable energy criterion is no longer required and 
the resulting Lyapunov function for the filter always has compact sub-level sets. In 
Section [5] we exploit this fact to prove geometric ergodicity of the signal-ensemble 
process for EnKF-AI and related methods. The ergodicity result is a type of stability 
which ensures that filter initialization errors are dissipated exponentially quickly by 
the method. 

In Section [7] we provide numerical evidence which supports the stability theory de¬ 
veloped for EnKF-AI and related methods. In particular, catastrophic filter divergence 
is completely eliminated by adaptive inflation, even when cheap unstable integrators 
are used in the forecast step. For such integrators, catastrophic Alter divergence is un¬ 
avoidable in EnKF and EnFK-CI. Moreover, adaptive inflation does not sacrifice any 
Altering skill, with EnKF-AI and EnKF-CAI displaying good accuracy when compared 
to benchmark accuracy scores, even in highly turbulent regimes. This suggests that 
adaptive inflation may in fact make filtering algorithms faster, since they allow for such 
cheap, typically unstable forecast integrators. 

In recent years, there have been many methods introduced featuring a type of 
adaptive inflation Altering. These methods are motivated by very different principles. 
For example, mm introduce a method that uses the norm of innovation to Alter the 
optimal parameter for multiplicative inflation. mm introduce two different ways 
to And the covariances of the system and observation noises, assuming they are not 
known. |15] designs a covariance relaxation method, so the posterior ensemble is pulled 
back to the forecast ensemble with a linear factor obtained from the innovation process. 
In [27] a modification of 3DVAR is proposed which relies on projection back to a stable 
region, this modification has been theoretically analyzed from the perspective of long¬ 
time accuracy. On one hand, these methods are quite different from the EnKF-AI 
method we introduced. With the exception of 1211, none have been rigorously proved 
to be stable in the sense of our theorems here. On the other hand, it is very possible 
that our theoretical framework here can be extended to these methods to build a 
stability framework, at least in some simple settings. The compelling reason here is 
that all these methods use the innovation sequence to decide the inflation, relaxation 
or noise covariance, which all grow with respect to O^- From this perspective, our 
framework here offers an explanation of the good performance of these data assimilation 
methods, as they are more resilient against Alter divergence; and when the Alter is on 
the edge of malfunctioning, the adaptive mechanisms can pull the ensemble back to 
more reasonable states. 


A Elementary claims 

Lemma A.l. By Young’s inequality, for any e > 0 and G the following holds: 

k + 2/P = -b |yp -b 2{x,y) < (1 -b -b (1 -b e~‘^)\y\^ . 

Lemma A.2. For any group of vectors vi,..., Vn, the following holds: 


arg mm 

V 




Vi -v = 


n 

-E 

n ^ 


Vi. 


k=l k=l 

Proof. One can verify this by using the KKT condition. 


□ 
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B Stability for ETKF-AI, EAKF-AI 


Proof for Theorem 4-5 We use the observed-unobserved notation introduced in Sec¬ 
tion 12.31 Notice that the mean update of (13.51) is the ensemble averaged version of 
(j3.ip without the artificial perturbations • Applying Jensen’s inequality 


_ 1 ^ 

En-l\Xn\^ < - 

k=l 

to (14.3p . we find the following with a constant Di, 

E„_i|X„|2 - 2pQ^\\HofEn-i\HUn\‘^ < Di. 

The corresponding part for ()4.4I1 will become 

|F„|2 < (1 + ^Ph)\Yn\^ + (1 + (I + HoC^H^)-\HoXn - 

Applying Jensen’s inequality to (14.51) . we find the second part above is bounded by 


+ HoC^H^)-\HoX^ - < K\\Hof 

So there is a constant D 2 such that 

En-l\Vn\^ - 2p^^\\HfEn-l\HUn\^ < (1 + + D^, 

where we use the fact that ||iJo|P = Then notice that the posterior covariance 

follows Cn Y Cn by (12.9p . Therefore using that 

K 

= a:|F„| 2+ tr(C'„) < JC|F„|2+ tr(CO, 

k=l 


and = K\Vri{^ +tr(C'n), we obtain from the previous inequality 

K K 

E„_1 ^ - 2Kp^^\\HfEn-l\HUn\^ < (1 + ^h)En-l ^ + D2. 

^=1 k=l 

Assumption 14.11 leads to 

K 1 ^ 1 

(1 + -Ph)En-i Y. < (1 - -Ph) Y + (1 + ifMKKh, 

k=l k=l 

and moreover 

K K 

E„_1 ^ yW|2-2i^po iF||X-i| Wnp < Y\vi%?+{l+-lih)KKh+D2. 

k=l k=\ 

Then notice a multiple of Assumption 14.11 is 

4A:/3-Vo < 4A:/3,:Vo (1 - h)\Un-i? + AKfi-^p^^\\HfKh. 
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The sum of the two previous two inequalities yields 


¥.n-l£n = + 4if/3-Vo ^(1 - 

k=l 

1 ^ 

< (1 - 2^^^) E(1 - mUn-1? + D, 

k=l 


with D := + (1 + ^f3h)KKh + D 2 . Then notice that (1 — ^ 

(1 — /3/i), therefore 

lEn-l^n < (1 ~ + D- 

□ 


C Ergodicity formulas 

Here we present the concrete formulas that are used in Section [5] for ETKF-AI and 
EAKF-AI. We omit unneccessary details concerning the transform and adjustment 
matrices, which can be found in [2l[3l[ll[7]. In both ESRF methods, the Markov kernel 
<I> : A X B{y) [0,1] with T := x R*^^^ and T := R'^ x x R"? is described by 

{Un-i, Eii\, ..., eS) ^ {Un, E^, ..., E^, Zn) . 

The deterministic step is given by the map r(17, E, Z) = {U, T^^),..., T^^^) where 

= V - CH^{I + HCH^)-\HV -Z) + (3.1) 

with E = v^k) and C = C = ^ - E) ® (E^^) - E) and 5^^) is 

the k-ih. column of the updated spread matrix S. S is given by ST{S) in the ETKF 
method, and A{S)S in the FAKE method, where S is the forecast spread matrix 

5 = (eW-E,...,eW-E). 

With the ETKF method, there are a few way to define the transformation matrix 
T(S'), one reasonable choice is taking the matrix square root 

T{S) = {Ik+{K- 1)~^S'^= {lK-{K-l)-^S'^H^{I+HCH^)-^HS'^f^‘^ . 

As for the EAKF method, the construction of the adjustment matrix A(S') is slightly 
more complicated as the following 

A{S) = QkG^{I + 

Here QAR is the SVD decomposition of S and DG is the diagonalization of {K — 
1)~^h?"Qhy, and f indicates pseudo inverse of a matrix. 

And for the ETKF-AI and EAKF-AI, everything is the same as their counterpart 
without adaptive inflation, except that we use the following in (13.11) 

~ 1 ^ - - 

^ = -K-^l ® or 

^ k=l 


32 




References 


[1] G. Evensen. The ensemble Kalman filter: Theoretical formulation and practical 
implementation. Ocean dynamics, 53(4):343-367, 2003. 

[2] C. H. Bishop, B. J. Etherton, and S. J. Majumdar. Adaptive sampling with the 
ensemble transform kalman filter, part i: Theoretical aspects. Mon. Weather Rev., 
129(3) :420-436, 2001. 

[3] J. L. Anderson. An ensemble adjustment Kalman filter for data assimilation. Mon. 
Weather Rev., 129(12):2884-2903, 2001. 

[4] A. J. Majda and J. Harlim. Filtering complex turbulent systems. Cambridge 
University Press, Cambridge, UK, 2012. 

[5] E. Kalnay. Atmospheric modeling, data assimilation, and predictability. Cam¬ 
bridge university press, 2003. 

[6] D. T. B. Kelly, K. J. Law, and A. M. Stuart. Well-posedness and accuracy of the 
ensemble Kalman filter in discrete and continuous time. Nonlinearity, 27:2579- 
2603, 2014. 

[7] X. T. Tong, A. J. Majda, and D. Kelly. Nonlinear stability and ergodicity of 
ensemble based kalman filters, submitted to Nonlinearity, 2015. 

[8] A.J. Majda and J. Harlim. Catastrophic filter divergence in filtering nonlinear 
dissipative systems. Comm. Math. Sci., 8:27-43, 2008. 

[9] G. Gottwald and A. J. Majda. A mechanism for catastrophic filter divergence 
in data assimilation for sparse observation networks. Nonlin. Processes Geophys, 
20:705-712, 2013. 

[10] D. Kelly, A. J. Majda, and X. T. Tong. Concrete ensemble kalman filters with 
rigorous catastrophic filter divergence. To apper in Proc. Natl. Acad. Sci. 

[11] J. L. Anderson. An adaptive covariance inflation error correction algorithms for 
ensemble filters. Tellus A, 59(210-224), 2007. 

[12] J L Anderson. Spatially and temporally varing adaptive covariance inflation for 
ensemble filters. Tellus A, 61A:72-83, 2009. 

[13] T. Berry and T. Sauer. Adaptive ensemble Kalman filtering of nonlinear systems. 
Tellus A, 65(20331), 2013. 

[14] Y. Zhen and J. Harlim. Adaptive error covariance estimation methods for ensemble 
Kalman filters. Journal of Computational Physics, 294:619-638, 2014. 

[15] Y. Ying and E. Zhang. An adaptive covariance relaxation method for ensemble 
data assimilation. Quarterly Journal of the Royal Meteorological Society, 2015. 

[16] J. Mandel, L. Cobb, and J. D. Beezley. On the convergence of the ensemble 
Kalman filter. Applications of Mathematics, 56(6):533-541, 2011. 

[17] E. Le Gland, Monbet V., and V. D. Tran. Large sample asymptotics for the 
ensemble kalman filter. Oxford Handbook of Nonlinear Filtering, 2010. 

[18] G. Evensen. Sampling strategies and square root analysis schemes for the EnKE. 
Ocean dynamics, 54(6):539-560, 2004. 


33 



[19] J. L. Anderson and S. L. Anderson. A Monte Carlo implementation of the non¬ 
linear filtering problem to produce ensemble assimilations and forecasts. Mon. 
Weather Rev., 127:2741-2758, 1999. 

[20] T. M. Hamill, C. Whitaker, and C. Snyder. Distance-dependent filtering of back¬ 
ground error covariance estimates in an ensemble Kalman filter. Mon. Weather 
Rev., 129:2776-2790, 2001. 

[21] R. Furrer and T. Bengtsson. Estimation of hign-dimensional prior and posterior 
covariance matrices in Kalman filter variants. Journal of Multivariate Analysis, 
98:227-255, 2007. 

[22] A. H. Jazwinski. Stochastic Processes and Filtering Theory. Academic Press, New 
York, 1972. 

[23] X. T. Tong and R. Van Handel. Ergodicity and stability of the conditional dis¬ 
tributions of nondegenerate markov chains. Ann. Appl. Probab., 22(4):1495-1540, 
2012. 

[24] S. Meyn and R. Tweedie. Markov chains and stochastic stability. Springer-Verlag, 
1993. 

[25] J. C. Mattingly and A. M. Stuart. Geometric ergodicity of some hypo-elliptic 
diffustion for particle motions. Markov Process. Related Fields, 8(2):199-214, 2002. 

[26] A. J. Majda, R. V. Abramov, and M. J. Grote. Information theory and stochastics 
for multiscale nonlinear systems, volume 25. American Mathematical Soc., 2005. 

[27] D. Sanz-Alonso and A. M. Stuart. Long-time asymptotics of the filtering distri¬ 
bution for partially observed chaotic dynamical systems. arXiv:1411.6510, 2014. 


34 



